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Abstract 



From any given Iterated Function System, a small set of balls that cover the 
fractal attractor can be simply determined. This gives a priori bounds on 
the region of space in which the attractor may be constructed. 



As Barnsley, Demko and others have shown [1, 2, 6, 4], one effective 
method for producing fractal shapes (in any number of dimensions) is with 
Iterated Function Systems (IFSs), using the “Chaos Game” algorithm (or 
some deterministic algorithm). This approach has been used for producing 
naturalistic shapes [4], finding fractal interpolants to given data [6, p. 274] 
and fractal approximations of given functions [8], and even for visualizing 
arbitrary discrete sequences [7]. Indeed, any contractive IFS will give an 
attractor (usually of fractal dimension); thus it is possible to generate IFSs at 
random to explore the graphical possibilities, as is done in some educational 
software [5]. Similarly, because the attractor depends continuously on the 
parameters in the IFS [1], a small data sets from any source could be encoded 
as IFSs for visualization. 

In implementing the IFS method, one important question is the prediction 
a priori of the region of space containing the fractal attractor. Without such 
a prediction, one could only approximately estimate the spatial extent based 
on calculating several points of the attractor, with no guarantee that these 
points are near the bounds. If as a result the portion of space represented 
in the computation of the attractor is too small, the result will not yield 
the whole attractor. If the portion of space represented is overly large, then 
much computational space is wasted, reducing the effective resolution of the 
computed attractor. Another concern is when the space has natural limits, 
for example the space of colors in Red-Green-Blue-space representable on a 
video monitor is limited (more or less) to a unit cube. Then the question 
becomes whether the attractor (or its projection onto the limited dimensions) 
will fit in the space. 

Here a simple algorithm is given to compute, directly from the IFS, a set 
of balls whose union contains the attractor as a subset (see Figure 1). The 
radii of the balls are minimal in a certain restricted sense. This gives reliable 
bounds on the region of space that must be considered in constructing the 
fractal. The method is general, independent of the particular space and 
metric. We first describe the set of balls, then show how to compute their 
radii and prove that the algorithm works, and lastly give a detailed example. 

An IFS consists of a set of n contraction mappings tc, : X -* X on 
a metric space X with metric d : X x X — ► TZ. (For the “Chaos Game” 
algorithm, probabilities pi are associated with each mapping; this idea has 
been extended to conditional probabilities [3]. Here only the former case is 
considered, where the attractor is independent of the non-zero p,.) Assume 



1 



that for each contraction w, the coutractivity ratio 0 < s, < 1 and the fixed 
point x, are known, where by definition s, satisfies d(w,(x), u»,(y)) < s,d(x, y) 
for all points x,y G X and x, satisfies x, = tc,(x t ). The action of the IFS IV' 
on a set S of points in X is defined as 

H’(S) = |J »,<S) (1) 

«=] 

where each contraction te, is applied to the set 5 in a pointwise sense. The 
attractor A is the set of points in X satisfying 

.4 = W(A) . (2) 

That is, the attractor consists of n smaller “copies" of itself. 

We seek to cover each of the n “copies” with a closed ball B, centered 
on the corresponding fixed point X;, so the radius r, must be chosen large 
enough that B, D u\(A). Call the union E of these balls the “envelope,” in 
that E D A by (2). Then relative to each x,, every point in the envelope will 
be within a distance R, = max_, ((/,_, + where = d(x,,x ; ), because for 
any point x in Bj , c/(x,, x) < d tJ + f/(x ; , x) < + rj. Applying w, to such a 

point x will give an image point y. where (/(x,, y) < s,c/(xi.x) < s,i?,. Hence 
if the radii r, are chosen to satisfy 



r, = + i'j ) (3) 

for i,j = 1 . . . n then B, will contain the image u’, (B) of the envelope and so 
E will contain its own image under the IFS: 

ED\V(E) . (4) 



Iterating the IFS from any starting set (E in particular) yields a sequence of 
sets that converges to the attractor. Since (4) implies E D W k (E) for any 
positive integer k, the envelope B, subject to (3), does indeed contain the 
attractor. But how can the r, be calculated from (3)? 

When n = 2 the radii can be determined algebraically. Solving the pair 
of equations (3) gives: 






r 2 



«l(l + * 2 ) j 

d ] 2 

1 — S1S2 

^2(1 + 5 l) J 

— di2 



( 5 ) 



2 



But for n > 2 there is apparently no closed-form general solution, and the r, 
must be found algorithmically. (If one wished to minimize calculating at the 
expense of overestimating, one could use r, = d maz s max /(l — s mar ), i.e., the 
envelope as if all fixed points were equidistant at the maximum separation 
and with all s, equal to the largest.) 

A natural approach for n > 2 is to start with the pairwise estimates 




Sj { 1 + Sj) 
1 - SiSj 



dtj 



max r„- 



( 6 ) 

(7) 



but in most cases the 7-j 1 * will not satisfy (3). The exceptional case is when 
r,j = rui for every i,j ^ k i, i.e., when for each ball all the pairwise 
estimates for that ball give the same size. (This case is not always apparent 
from the attractor: Figure 1 shows such an example.) Otherwise, some of 
the rj 1 * will be too small to contain some images ir,-(i?j 1) ) of the other balls. 
Then the obvious iterative scheme to try is 

rj* +1) = max «<(</„ + rf ] ) , t, j = 1 . . . n . (8) 



Because this approach never overestimates the radii (rj A> < r,) and the iter- 
ates are nondecreasing, the algorithm must converge. What is not so obvious 
is that this process always suceeds in at most n — 1 iterations, as shown be- 
low. (This would not be true without (7).) In fact, there is a direct algorithm 
(not iterative) that is more efficient. 

The key idea that allows the direct algorithm is that the distances d,j can 
be rescaled to account for the contraetivities s,, and the scaled distances 
can be used to order the contractions w,. Let 



, (1 + 5,’)(1 + Sj) 

,j ~ 1 - SiSj 



i,j = 1 . . .n 



(9) 



(While Dij is clearly symmetric and non-negative, it is not a metric because 
it doesn’t satisfy the triangle inequality.) Now reorder (and relabel) the w, 
by decreasing maximum scaled distance, so that 



i < j => max > max Dji , i,j, kj = 1 . . . n (10) 
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In the new order, use the pairwise formula (5) for the first two radii. Then 
proceed in order based on the previous results, letting 

r, = s, max(d,j + rj) , i = 3 . . . n . (11) 

This is the direct algorithm, which, as shown below, solves (3); an imple- 
mentation in the C programming language is given in the Appendix. 

First note that in the exceptional case mentioned above (after 7), all the 
Djj (for i ^ j) are equal. In this case, the direct algorithm will obtain the 
correct r, regardless of the order in which they are computed. If not all the 
Djj are equal, then some of the r, will need to be larger than the pairwise 
estimates rj 1 *, and hence larger than the r tJ in (6), so in general (3) implies 

( 12 ) 

for i,j = 1 . . . n. 

For the general case, the proof is by induction, showing that each new r, 
computed requires no adjustment of those previously computed. Clearly rj, 
r -2 from (5) satisfy (3) for the subset i,j = 1,2. Now assume the first m — 1 
radii, in the order (10). satisfy (3), and hence (12), for i,j = l...m — 1. 
Choose r„, by (11). and let k be the value of the index j in (11) for which 
the maximum is achieved. Then by (12) and the ordering (10) 

r* > rr-^b, • (13) 

1 + s k 

Algebraic manipulation of (13) gives 

? k ^ *t[(l d" $m)dkm d" SpiU) = Sk((/k m d" 1 m ) (Id) 



so the new r m requires no alteration of r k . 
Similarly, for ? ^ A\ i < m 



r, > 

> 



J~ D 

d- S 



s, ( 1 •+■ s m ) 

1 ^ vi 



since s, < 1. Combining (15) with 



r, > Sj{d ik + r k ) 



(15) 

(16) 
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yields 



ri ^ ^t[(I d” ^m)djni "(■ -(- r^.)] 

^ S,[d trn 4* Sm(d);m d" rit)] = ^i(d|'rn 4- 7%,, ) (10 

by the triangle inequality. Hence the new r m requires no adjustment of any r, 
for i < m, and (3) is satisfied for i,j = 1 . . . m; this completes the proof. Note 
also that because of how the direct algorithm works, the iterative algorithm 
will compute at least two of the r,- (/q , r 2 in the order (10)) in the initial step, 
and will find at leqst one of the other r, at each successive step, and so can 
take at most n — 1 iterations to arrive at the answer. 

Can radii smaller than these r, be used and still have the B, cover the 
attractor? For any particular IFS. the answer is probably yes (as illustrated 
in Figure 2). The approach given here uses only minimal information about 
the IFS: the ordering of the based on the maximal scaled distances; the 
s,-; and for each i one determining distance (where the maximum in (11) 
is achieved). Using more information it may be possible to reduce the size 
of the B{. But if one considers the set of all the IFSs for which the direct 
algorithm yields the same r, in the same way (i.e., same ordered s, and same 
n — 1 determining distances), then the r, are minimal for that set of attractors 
(see Figure 3). In fact, one can construct one member of that set such that 
each image u’,(.4) of the attractor includes a point at a distance r, from the 
fixed point x,. 

To construct this IFS. let X = 1Z with the Euclidean metric d(x,y) = 
|.r — y\. Let u'i(x) = — s,x 4- (1 + s,).t,. and let .Tj = 0, x 2 = d\o, say. Then 
the attractor includes the extremal points x e \ = .rj — jq, x e 2 = x 2 + r 2 . since 
x e i = ti’i (x c2 ) and x c 2 = w 2 (x c i ). Place each succedingx TO at the determining 
distance d m k from the determining point x*. in the opposite direction from 
X(h . (Figure 4 illustrates the construction.) Then the attractor will include 
x em = x„, ± = iv m (x e k)- Thus for this one-dimensional attractor A, each 

image Wi(A) will include a point (x et ) a distance r, from x,-, so no smaller r, 
would suffice. 

As an example, consider the now familiar black spleenwort fern fractal of 
[4]. The IFS for the fern (in two dimensions) consists of affine contractions, 
each of which has the form 




o 



or more compactly 



w,(x) = M,x + b, . (19) 

where M, is the matrix and &,• is the offset vector. The various constants are 
given in [ 4 ], but in terms of scaling and rotating each axis, using p, q, 9, <\> , 
where a = p cos 9, b = —q sin0, c = psin0, d = qcos<t>. The following table 
is adapted from [ 4 , p. 1977]: 



Map 


Scalings 


Rotations 


Translations 


i 


Pi 


9. 


9i <p, 


e, 


fi 


1 


0 


0.16 


0 


0 


0 


0 


2 


0.85 


0.85 


-2.5 


-2.5 


0 


1.6 


3 


0.3 


0.34 


49 


49 


0 


1.6 


4 


0.3 


0.37 


120 


-50 


0 


0.44 



where angles are given in degrees. The contractivity ratio s, for an affine 
map Wi is the largest singular value of the matrix Mi. In the first three maps 
above, both axes rotate together, and so s,- is the larger of p,, <p. In w 4 . the 
differential rotation causes a skewing effect, and the singular values of M 4 
must be found. The simplest way for a real 2x2 matrix is first ti> factor out 
a pure rotation to give a symmetric matrix (5), then diagonalize it to find 
its eigenvectors (Ai. A 2 ) as shown below: 



S = 




= arctan 



' c-fr N 

a + d . 



3 = 



[ ^ 


0 1 


_ / 


[ 0 


A* j 


~ \ 



/ COSO 


si no 


\ —sin a 


COS Q 


1 / 

- arctan 
2 \ 


' 2 h ' 


.9~ k . 


( cos 3 


sin 3 


1 — sin;? 


cos 3 



M 




( 20 ) 

( 21 ) 

( 22 ) 

(23) 



Then s — max(|Ai|, | A 2 1 ) . This approach also has a nice geometrical inter- 
pretation: the effect of multiplying a vector x by M is to take components 
of x in the eigenvector directions, which are at an angle /? relative to the 
coordinate axes, scale each component by the corresponding A, and rotate 
the resulting vector by a. 

Proceeding as above gives s 4 = 0.379. The fixed point x, for each map 
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can be found by solving 



(I - Mi) x, = b, (24) 

This gives the necessary starting information, summarized as “Input” in the 
table below. For affine maps in higher dimensions, the contractivity ratios 
are found by singular value decomposition, but for nonlinear maps the ratios 
and fixed points may be more difficult to find. 

Running the direct algorithm program (“envelope”) from the Appendix 
on this input gives the following results: 

7, envelope 4 

Map 1. Enter scale, x, y: .16 0 0 

Map 1: s = 0.160000, x = 0.000000, y = 0.000000 
Map 2. Enter scale, x, y: .85 2.45967 10.004734 

Map 2: s = 0.850000, x = 2.459670, y = 10.004734 
Map 3. Enter scale, x, y: .34 -0.601889 1.883961 

Map 3: s = 0.340000, x = -0.601889, y = 1.883961 
Map 4. Enter scale, x, y: 0.379216 0.155336 0.630251 

Map 4: s = 0.379216, x = 0.155336, y = 0.630251 
radii in sorted order [orig order] (sorted link): 
rl[2](->2): 16.700212 
r2[4](->l): 9.993765 
r3 [3] (->1) : 8.628835 
r4 [1] (->1) : 4.320459 
7 . 



These results are illustrated in Figure 5, and tabulated below (including the 
D values used in re-ordering and the determining distances): 



Input 


Output 


i 


Si 


•T. 


y, 


D- order 


Dmax 


ddet 


ri 


i 


0.16 


0 


0 


4 


D\ 2 ' 25.59 


d 12 


10.30 


4.32 


2 


0.85 


2.460 


10.005 


1 


Z?24* 36.35 


d 24 


9.65 


16.70 


3 


0.34 


-0.602 


1.884 


3 


D32' 30.26 


d 32 


8.68 


8.63 


4 


0.379 


0.155 


0.630 


2 


D42' 36.35 


d 42 


9.65 


9.99 



the same pair that gives the maximum D, that is not always the case.) Then, 
If we had no idea how big the fern attractor was, we could use a computational 
space extending from x mi „ = .r 2 — r 2 = —14.24 to x max = x 2 + r 2 = 19.16 



and y mtn = y 4 — r 4 = —9.36 to y maz = y 2 + r 2 = 26.70 to contain the entire 
envelope. As it turns out, this is far more space than necessary for the fern 
itself, but there are many other IFSs, equivalent as far as the direct algorithm 
is concerned, with much larger attractors (e.g., what if 0 2 = = 177.5 

instead). 

To summarize, given any IFS (along with the contractivities and fixed 
points of each of its constituent contraction mappings), an envelope can be 
constructed of one ball for each map, centered on the corresponding fixed 
point. (In the case of affine maps in two dimensions, an explicit procedure 
for finding the contractivities and fixed points was given.) We have proven 
that the radii of the balls can be calculated by a simple algorithm (direct 
or iterative) such that the envelope covers the entire attractor. The spatial 
extent of the envelope thus gives a reliable bound on that of the attractor. (In 
addition, if the balls are disjoint, the attractor is totally disconnected.) While 
the radii found by the direct algorithm may not be minimal for the particular 
IFS. they are minimal for the set of all IFSs with equivalent information (in 
the sense described above). 
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APPENDIX 

/* 

The following program in C implements the direct algorithm 
for determining the "envelope" of an attractor of an Iterated 
Function System on R“2, given the contractivities si and the 
fixed points (xi,yi). 

written by David Canright, March 1993. 

*/ 

#include <stdio.h> 
finclude <math.h> 

int npts, i, j, m, n, index [64], link [64]; 
double d[64] [64] , x[64] , y[64], s[64], r[64], Dmax[64], 
t, tmax, dx, dy; 

int input (int argc, char *argv[]); 
main(int argc, char *argv[]) { 

npts = input (argc, argv); /* get si, xi, yi */ 

/* compute distances dij & maximal scaled Dij */ 
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for (i = 1; i <= npts; i++) Dmax[i] =0.; 

for (i = 1; i <= npts; i++) { 

for (j = i+1; j <= npts; j++) { 

dx = x[i]-x[j] ; dy = y[i]-y[j]; 
d[i] [j] = d[j][i] = t = sqrt(dx*dx+dy*dy) ; 
t = ( 1 . +s [i] ) * ( 1 . +s [ j] ) / ( 1 . -s [i] *s [j ] ) * t; 
if (t > Dmax[i]) Dmax[i] = t; 

if (t > Dmax[j]) Dmax[j] = t; 

> 

> 

/* Sort by scaled distances; index points to old order */ 
index [1] = 1; 

for (i = 2; i <= npts; i++) { 

for (m = i; m > 1 kk Dmax[i] > Dmax [index [m-1] ] ; m — ) 
index[m] = index[m-l]; 
index [m] = i ; 

> 

/* Direct algorithm; link points to determining distance */ 
i = index[l]; j = index[2]; link[l] = 2; link[2] = 1; 
r [1] = ( s[i]/(l .+s[i] ) ) * Dmax[i]; 
r[2] = ( s [j ] / ( 1 .+s[j] ) ) * Dmax[j]; 
for (m=3; m <= npts; m++) { 
i = index [m] ; 
tmax = 0. ; 

for (n = 1; n < m; n++) { 
j = index [n] ; 

if ( (t = d[i][j] + r[n]) > tmax ) 

{ tmax = t; link[m] = n; > 

> 

r[m] = s[i] * tmax; 

> 

printf ("radii in sorted order [orig order] (sorted link):\n"); 
for (i = 1; i <= npts; i++) 

printf (" r*/.d[*/.d] (->’/.d) : 7, f \n" , i , index [i] , link [i] ,r [i] ); 
return(O) ; 

> 
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/* Input function: gets si, xi, yi, or uses random numbers */ 

/* optional arguments: number of maps, seed for random */ 
int input (int argc, char *argv[]) { 
double norm; 

char line [81], getmore = 1; 

npts =3; /* default */ 

if (argc > 1) sscanf (argv[l] ,'"/d",ftnpts); 
if (argc > 2) sscanf (argv[2] , "*/,d" ,fti) ; 
srand(i); 

norm = l./(MAXINT); /* machine-dep. const., to normalize */ 
for (i = 1; i <= npts; i++) { 

s[i] = norm*rand(); /* random by default */ 
x[i] = norm*rand(); 
y[i] = norm*rand(); 

/* get numbers from stdin until blank line, then random */ 
if (getmore) { 

printfC'Map */,d. Enter scale, x, y: ",i); 
gets(line) ; 
if (line[0]) { 

sscanf (line , '"/.lf'/lf’/.lf " ,s+i ,x+i ,y+i) ; 

if (s[i]<0. )s[i] = -s[i]; /* enforce 0 <= s < 1 */ 

while(s [i] >=1 . ) s[i] *= 0.1; 

} 

else getmore = 0; 

> 

printfC'Map */,d: s = 7,f, x = */,f, y = */f \n" , i,s [i] ,x[i] ,y [i] ) ; 

> 

return(npts) ; 

> 
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FIGURE CAPTIONS 



Figure 1. An attractor of an IFS is shown with its envelope of three disks, as 
computed by the direct algorithm. (This IFS uses affine maps, with sj = |, 
«2 = 5 , 53 = §, xi = (0,0), x 2 = (4,0), and x 3 = (0,3). 

Figure 2. For an equilateral Sierpinski’s Triangle (where w t -(x) = |x+ ^x,), 
the r, = dij by this method; in this particular case the radii could be half as 
large. 

Figure 3. Same s, and x, as the previous figure, but here Wi(x) = — ^x+ |x,; 
in this case the r, found above are minimal. 

Figure 4. A one-dimensional attractor constructed from the following or- 
dered data: si = 5 , s 2 = s 3 = ;• d 12 = 1, and The r, are minimal 

for such attractors. 

Figure 5. Barnsley’s fern [4] and its envelope (see text). 
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